本日は地震観測データとその他の観測データを組合せて可視化する手法を学びます.
地震観測データには,アメリカ地質調査所(United States Geological Survey; USGS) Earthquake Hazards Program(https://earthquake.usgs.gov/earthquakes/)が公開するデータを用います.
USGSの地震観測データと組み合わせるその他の観測データには, テキサス大学地球物理研究所PLATESプロジェクト(http://ig.utexas.edu/marine-and-tectonics/plates-project/)が公開するプレート境界線データとスミソニアン博物館国立自然史博物館GLOBAL VOLCANISM PROGRAM(http://www.volcano.si.edu/)が公開する火山データを用います.
本日の演習では,緯度(latitude)・経度(longitude)という情報を基に複数のデータを組み合わせる手法の習得を目標とします.
まず,テキサス大学地球物理研究所PLATESプロジェクト(http://ig.utexas.edu/marine-and-tectonics/plates-project/)が公開するプレート境界線データを可視化します。
プレート境界線とはプレートテクトニクス(プレート理論)により説明されるプレートとプレートの境目のことです.詳しくは以下リソースを参照してください.
http://ja.wikipedia.org/wiki/プレートテクトニクス
テキサス大学地球物理研究所PLATESプロジェクトのデータ公開ページ(http://www-udc.ig.utexas.edu/external/plates/data.htm)からプレート境界線データをダウンロードします.
ページの中頃に「PLATES_PlateBoundary_ArcGIS.zip」というファイルがあるので,それをダウンロードします。
Rの作業ディレクトリにディレクトリ「plate」を作成し, そのディレクトリに「PLATES_PlateBoundary_ArcGIS.zip」を展開してください.
Rで上記処理をさせる場合は以下となります.
# fileがなければダウンロードする
if(!file.exists("PLATES_PlateBoundary_ArcGIS.zip")){
download.file(url="http://www-udc.ig.utexas.edu/external/plates/data/plate_boundaries/PLATES_PlateBoundary_ArcGIS.zip",destfile = "PLATES_PlateBoundary_ArcGIS.zip")
}
# directoryがなければファイルを展開する
if(!dir.exists("plate")){
unzip("PLATES_PlateBoundary_ArcGIS.zip",exdir = "plate")
}ダウンロードしたプレート境界線データはESRI Shapeフォーマットで格納されています. Shape形式はGIS(地理情報システム)で比較的よく利用されているファイルフォーマットです. 詳しくは Esri Japanのシェープファイルについて等を参照してください.
Shapeファイルを処理するため maptools ライブラリをインストールします. また,データフレーム変換用のライブラリ broomもインストールします.
#ライブラリの読み込み.なければインストールする.
if(!require(ggplot2)){
install.packages("ggplot2")
library(ggplot2)
}
if(!require(maps)){
install.packages("maps")
library(maps)
}
if(!require(maptools)){
install.packages("maptools")
library(maptools)
}
if(!require(broom)){
install.packages("broom")
library(broom)
}
#游ゴシック体を使う
if(.Platform$OS.type=="windows")
windowsFonts(yugo=windowsFont("Yu Gothic"))
if(capabilities("aqua"))
quartzFonts(yugo=quartzFont(rep("YuGo-Medium",4)))それでは早速Shapeファイルを描画しましょう.
展開したShapeファイルを読み込み,ggplotで発散型境界を描画します.
ggplot2 は data.frame でないと処理できないので, maptoolsの処理結果(SpatialLinesDataFrame)をそのまま利用することができません. そこで,broom パッケージのtidy 関数を利用して SpatialLinesDataFrame を data.frame に変換ます.
それでは早速下記コードを実行してください.
# shapeファイルの読み込み
ridge <-readShapeSpatial("plate/ridge.shp")
# オブジェクトのクラスを確認する
class(ridge)## [1] "SpatialLinesDataFrame"
## attr(,"package")
## [1] "sp"
ridge.fort<-broom::tidy(ridge)
class(ridge.fort)## [1] "data.frame"
g <- ggplot() + geom_path(data= ridge.fort, aes(x=long, y=lat, group=group),color="red")+theme_bw(base_family = "yugo")+ggtitle("発散型境界 (data by UTIG Plate Proj.)")
g発散型境界線が描画されました.
世界地図と合わせて描画してみましょう.
g+borders("world")つぎにトランスフォーム型境界を描画します。
transform <-readShapeSpatial("plate/transform.shp")
transform.fort<-broom::tidy(transform)
g <- ggplot() + geom_path(data= transform.fort, aes(x=long, y=lat, group=group),color="red")+theme_bw(base_family = "yugo")+ggtitle("トランスフォーム型境界 (data by UTIG Plate Proj.)")
gこれも世界地図と合わせて描画してみましょう.
g+borders("world")つぎに収束型境界を描画します。
trench <-readShapeSpatial("plate/trench.shp")
trench.fort<-broom::tidy(trench)
g <- ggplot() + geom_path(data= trench.fort, aes(x=long, y=lat, group=group),color="red")+theme_bw(base_family = "yugo")+ggtitle("収束型境界 (data by UTIG Plate Proj.)")
g収束型境界についても,世界地図と合わせて描画してみましょう.
g+borders("world")それでは発散型,トランスフォーム型,収束型境界線を合わせて描画します.
ridge <-readShapeSpatial("plate/ridge.shp")
ridge.fort<-broom::tidy(ridge)
transform <-readShapeSpatial("plate/transform.shp")
transform.fort<-broom::tidy(transform)
trench <-readShapeSpatial("plate/trench.shp")
trench.fort<-broom::tidy(trench)
g <- ggplot()
g <- g + theme_bw(base_family = "yugo")+ggtitle("プレート境界 (data by UTIG Plate Proj.)")
#発散型境界の描画
g <- g + geom_path(data= ridge.fort, aes(x=long, y=lat, group=group),color="red")
#トランスフォーム型境界の描画
g <- g + geom_path(data= transform.fort, aes(x=long, y=lat, group=group),color="blue")
#収束型境界の描画
g <- g + geom_path(data= trench.fort, aes(x=long, y=lat, group=group),color="green")
g世界地図と合わせます.
g + borders(database = "world")地震はプレート境界線で生じる歪が原因となり発生するといわれています.
日本で地震が多い理由も納得できますね.
アメリカ地質調査所(United States Geological Survey; USGS) のEarthquake Hazards Program(https://earthquake.usgs.gov/earthquakes/)の地震観測データとテキサス大学地球物理研究所PLATESプロジェクトのプレート境界データを合成します.
まずは地震観測データを描画します. データは前回の可視化で用いた「20110311JSTquery.csv」を使います.
csv<-read.csv("20110311JSTquery.csv",stringsAsFactors = FALSE)
#地震をプロット(色で深度,サイズでマグニチュードを表現)
g<-ggplot(csv,aes(x=longitude,y=latitude,color=depth,size=mag))
g<-g+geom_point()
g<-g+borders("world")
g<-g+theme_bw(base_family = "yugo")
g<-g+ggtitle("地震データ(data by USGS Earthquake Hazards Prog. Earthquake Hazards Prog.)")
g<-g+ylab("緯度")+xlab("経度")
#深度が浅ければ暗く,深ければ明るくする
g<-g+scale_color_continuous(low=rgb(1,0,0),high=rgb(1,1,0))
#点のサイズを0から15の範囲とする
g<-g+scale_size_continuous(range=c(0,6))
gここまでは前回と同じですね.
それでは地震観測データとプレート境界線データと合わせて可視化してみましょう.
ポイントはこれまで描画データを指定していたggplot()関数では何も指定せずに, 点や線を描画する geom_point や geom_path で描画データを指定している点です. こうすることで複数のデータを用いた可視化が可能となります(ggplot 関数でデータを指定すると全体でデータが共有されることになる).
# 地震観測データの読み込み
csv<-read.csv("20110311JSTquery.csv",stringsAsFactors = FALSE)
# プレート境界線データの読み込み
ridge <-readShapeSpatial("plate/ridge.shp")
ridge.fort<-broom::tidy(ridge)
transform <-readShapeSpatial("plate/transform.shp")
transform.fort<-broom::tidy(transform)
trench <-readShapeSpatial("plate/trench.shp")
trench.fort<-broom::tidy(trench)
# 描画データを指定しない
g<-ggplot()
# geom_point で描画するデータ(地震観測データ)を指定している.
g<-g+geom_point(data = csv,mapping = aes(x=longitude,y=latitude,color=depth,size=mag))
# 世界地図の描画
g<-g+borders("world")
g<-g+theme_bw(base_family = "yugo")
g<-g+ggtitle("地震データ(data by USGS Earthquake Hazards Prog.)")
g<-g+ylab("緯度")+xlab("経度")
#深度が浅ければ暗く,深ければ明るくする
g<-g+scale_color_continuous(low=rgb(1,0,0),high=rgb(1,1,0))
#点のサイズを0から15の範囲とする
g<-g+scale_size_continuous(range=c(0,6))
# タイトルを上書き
g <- g +ggtitle("プレート境界と地震データ(data by UTIG Plate Proj. and USGS)")
#発散型境界の描画(データを指定している点に注目)
g <- g + geom_path(data= ridge.fort, aes(x=long, y=lat, group=group),color="red")
#トランスフォーム型境界の描画(データを指定している点に注目)
g <- g + geom_path(data= transform.fort, aes(x=long, y=lat, group=group),color="blue")
#収束型境界の描画(データを指定している点に注目)
g <- g + geom_path(data= trench.fort, aes(x=long, y=lat, group=group),color="green")
gプレート境界線近辺で地震が発生していることがわかります.
地震発生メカニズム(気象庁,地震発生のしくみ)を実際のデータで確認することができました.
つぎにプレート境界線データと火山データを組合せて可視化してみます.
火山もプレートテクトニクスに関連が強い地形です.(参考:Wikipdia,火山)
火山データには,スミソニアン博物館国立自然史博物館GLOBAL VOLCANISM PROGRAM(http://www.volcano.si.edu/)が公開する火山データを用います.
火山データのダウンロードするため,スミソニアン博物館国立自然史博物館GLOBAL VOLCANISM PROGRAMサイト(http://www.volcano.si.edu/)へアクセスします.
「Database」 > 「HOLOCENE SPREADSHEET」をクリックします.
「GVP_Volcano_List.xls」というファイル名でダウンロードが始まります.
ダウンロードが完了したら,「GVP_Volcano_List.xls」をCSVフォーマットに変換します.
フォーマット変換にはExcelを利用します.Excelでファイルを開き,「名前を付けて保存」でCSVフォーマットとして保存します.
変換後のCSVファイル「GVP_Volcano_List.csv」をRの作業ディレクトリに移動します.
それではさっそく火山データを読み込んでみましょう.
#データを読み込む
volc <- read.csv("GVP_Volcano_List.csv",head=TRUE,skip = 1,stringsAsFactors = FALSE)
#構造を確認する
str(volc)## 'data.frame': 1520 obs. of 13 variables:
## $ Volcano.Number : int 210010 210020 210030 210040 211001 211003 211004 211010 211020 211030 ...
## $ Volcano.Name : chr "West Eifel Volcanic Field" "Chaine des Puys" "Olot Volcanic Field" "Calatrava Volcanic Field" ...
## $ Country : chr "Germany" "France" "Spain" "Spain" ...
## $ Primary.Volcano.Type: chr "Maar(s)" "Lava dome(s)" "Pyroclastic cone(s)" "Pyroclastic cone(s)" ...
## $ Activity.Evidence : chr "Eruption Dated" "Eruption Dated" "Evidence Credible" "Eruption Dated" ...
## $ Last.Known.Eruption : chr "8300 BCE" "4040 BCE" "Unknown" "3600 BCE" ...
## $ Region : chr "Mediterranean and Western Asia" "Mediterranean and Western Asia" "Mediterranean and Western Asia" "Mediterranean and Western Asia" ...
## $ Subregion : chr "Western Europe" "Western Europe" "Western Europe" "Western Europe" ...
## $ Latitude : num 50.2 45.8 42.2 38.9 43.2 ...
## $ Longitude : num 6.85 2.97 2.53 -4.02 10.87 ...
## $ Elevation..m. : int 600 1464 893 1117 500 800 949 458 1281 789 ...
## $ Dominant.Rock.Type : chr "Foidite" "Basalt / Picro-Basalt" "Trachybasalt / Tephrite Basanite" "Basalt / Picro-Basalt" ...
## $ Tectonic.Setting : chr "Rift zone / Continental crust (>25 km)" "Rift zone / Continental crust (>25 km)" "Intraplate / Continental crust (>25 km)" "Intraplate / Continental crust (>25 km)" ...
#要約を確認する
summary(volc)## Volcano.Number Volcano.Name Country
## Min. :210010 Length:1520 Length:1520
## 1st Qu.:261148 Class :character Class :character
## Median :300020 Mode :character Mode :character
## Mean :296768
## 3rd Qu.:342123
## Max. :600000
## Primary.Volcano.Type Activity.Evidence Last.Known.Eruption
## Length:1520 Length:1520 Length:1520
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
## Region Subregion Latitude Longitude
## Length:1520 Length:1520 Min. :-78.50 Min. :-179.97
## Class :character Class :character 1st Qu.: -7.13 1st Qu.: -78.15
## Mode :character Mode :character Median : 13.84 Median : 38.39
## Mean : 14.05 Mean : 23.85
## 3rd Qu.: 42.00 3rd Qu.: 139.04
## Max. : 85.61 Max. : 179.58
## Elevation..m. Dominant.Rock.Type Tectonic.Setting
## Min. :-4200.0 Length:1520 Length:1520
## 1st Qu.: 680.8 Class :character Class :character
## Median : 1462.0 Mode :character Mode :character
## Mean : 1677.0
## 3rd Qu.: 2344.8
## Max. : 6879.0
#データを確認する
head(volc)## Volcano.Number Volcano.Name Country Primary.Volcano.Type
## 1 210010 West Eifel Volcanic Field Germany Maar(s)
## 2 210020 Chaine des Puys France Lava dome(s)
## 3 210030 Olot Volcanic Field Spain Pyroclastic cone(s)
## 4 210040 Calatrava Volcanic Field Spain Pyroclastic cone(s)
## 5 211001 Larderello Italy Explosion crater(s)
## 6 211003 Vulsini Italy Caldera
## Activity.Evidence Last.Known.Eruption Region
## 1 Eruption Dated 8300 BCE Mediterranean and Western Asia
## 2 Eruption Dated 4040 BCE Mediterranean and Western Asia
## 3 Evidence Credible Unknown Mediterranean and Western Asia
## 4 Eruption Dated 3600 BCE Mediterranean and Western Asia
## 5 Eruption Observed 1282 CE Mediterranean and Western Asia
## 6 Eruption Observed 104 BCE Mediterranean and Western Asia
## Subregion Latitude Longitude Elevation..m.
## 1 Western Europe 50.170 6.85 600
## 2 Western Europe 45.775 2.97 1464
## 3 Western Europe 42.170 2.53 893
## 4 Western Europe 38.870 -4.02 1117
## 5 Italy 43.250 10.87 500
## 6 Italy 42.600 11.93 800
## Dominant.Rock.Type
## 1 Foidite
## 2 Basalt / Picro-Basalt
## 3 Trachybasalt / Tephrite Basanite
## 4 Basalt / Picro-Basalt
## 5 No Data (checked)
## 6 Trachyte / Trachydacite
## Tectonic.Setting
## 1 Rift zone / Continental crust (>25 km)
## 2 Rift zone / Continental crust (>25 km)
## 3 Intraplate / Continental crust (>25 km)
## 4 Intraplate / Continental crust (>25 km)
## 5 Subduction zone / Continental crust (>25 km)
## 6 Subduction zone / Continental crust (>25 km)
#富士山のデータを確認する
volc[volc$Volcano.Name=="Fujisan",]## Volcano.Number Volcano.Name Country Primary.Volcano.Type
## 598 283030 Fujisan Japan Stratovolcano
## Activity.Evidence Last.Known.Eruption Region
## 598 Eruption Observed 1708 CE Japan, Taiwan, Marianas
## Subregion Latitude Longitude Elevation..m. Dominant.Rock.Type
## 598 Honshu 35.361 138.728 3776 Basalt / Picro-Basalt
## Tectonic.Setting
## 598 Subduction zone / Continental crust (>25 km)
火山データの位置が緯度(Latitude),経度(Longitude)に記録されていることがわかります.
それでは火山データを可視化してみましょう.ここでは散布図 geom_point を用いていますが,単なる点だと火山ぽくないので,「shape=24」を指定することで点を三角に変更しています.
#火山をプロットする shape=24で点のサイズを変更する
ggplot()+geom_point(data=volc,mapping = aes(x=Longitude,y=Latitude),shape=24,color="#FF0000",fill="#FF0000")+theme_bw(base_family = "yugo")+ggtitle("火山データ\n(data by Smithsonian NMNH Global Volcanism Program,)")geom_point で指定可能な shape には他にも以下のような形があります.
g<-ggplot(data=data.frame(x=0:24))+scale_shape_identity()
g<-g+geom_point(mapping=aes(x=x,y=0,shape=x),size=5)
g<-g+geom_point(mapping=aes(x=x,y=1,shape=x),size=5,fill="red")
g<-g+xlab("shape types")
g今回はshapeの24番(色塗り可能な三角形)を指定することで火山を表現しています.
スミソニアン博物館国立自然史博物館GLOBAL VOLCANISM PROGRAMの火山データにはその他にも様々な情報が記録されています.例えば,「Elevation..m.」には火山の標高が記録されています.
可視化してみると次のようになります.
volc <- read.csv("GVP_Volcano_List.csv",head=TRUE,skip = 1)
#火山をプロットする shape=24で点のサイズを変更する
ggplot()+geom_point(data=volc,mapping = aes(x=Longitude,y=Latitude,color=Elevation..m.,fill=Elevation..m.),shape=24)+theme_bw(base_family = "yugo")+ggtitle("火山データ\n(data by Smithsonian NMNH Global Volcanism Program)")それでは火山データとプレート境界線データを合わせて描画してみましょう.
まず,火山データと世界地図を描画します.
volc <- read.csv("GVP_Volcano_List.csv",head=TRUE,skip = 1)
#火山をプロットする shape=24で点のサイズを変更する
g<-ggplot()
g<-g+borders(database = "world")
g<-g+geom_point(data=volc,mapping = aes(x=Longitude,y=Latitude),shape=24,color="#FF0000",fill="#FF0000")
g<-g+theme_bw(base_family = "yugo")
g<-g+ggtitle("火山データ\n(data by Smithsonian NMNH Global Volcanism Program,)")
gこの描画結果にプレート境界線を追加すればよいわけです.
地震観測データとプレート境界線データの描画を参考にすると, 下記のようにコードを書くことができると思います.要は描画単位(geom_pointやgeom_pathごと)に可視化するデータを指定すればよいのです.
#プレート境界線データの読み込み
ridge <-readShapeSpatial("plate/ridge.shp")
ridge.fort<-broom::tidy(ridge)
transform <-readShapeSpatial("plate/transform.shp")
transform.fort<-broom::tidy(transform)
trench <-readShapeSpatial("plate/trench.shp")
trench.fort<-broom::tidy(trench)
#火山データの読み込み
volc <- read.csv("GVP_Volcano_List.csv",head=TRUE,skip = 1)
g<-ggplot()
#地図の描画
g<-g+borders(database = "world")
#発散型境界の描画
g <- g + geom_path(data= ridge.fort, aes(x=long, y=lat, group=group),color="red")
#トランスフォーム型境界の描画
g <- g + geom_path(data= transform.fort, aes(x=long, y=lat, group=group),color="blue")
#収束型境界の描画
g <- g + geom_path(data= trench.fort, aes(x=long, y=lat, group=group),color="green")
#火山の描画
g<-g+geom_point(data=volc,mapping = aes(x=Longitude,y=Latitude),shape=24,color="#FF0000",fill="#FF0000")
g<-g+theme_bw(base_family = "yugo")
g<-g+ggtitle("火山データ\n(data by Smithsonian NMNH Global Volcanism Program,)")
gこのように火山もプレート境界線近辺に多くできていることがわかります.
本日は地球観測データと科学データを複数組合せて可視化する手法を学びました. 具体的には,緯度(latitude)・経度(longitude)をキーに,地震観測データ,プレート境界線データ,火山データを関連付けて可視化しました.このように複数のデータを組合せて可視化するにはデータを関連付けるためのキーが必要となります.今回の演習で感じてもらえたと思いますが,緯度・経度は複数データを組み合わせる際のキー候補としやすい情報の一つです.
次回は今回の内容に加えて地理空間情報の可視化に取り組みます.
以下の問題を解いてください. その際に用いたRスクリプトファイル,結果として得られた画像ファイルをZIPでまとめて提出してください.
※Rスクリプトファイルは可視化毎に作成してください.可視化に必要な内容のみを含めるようにしてください. ※データファイルは不要です.
ヒント:
可視化例:
ヒント:
例:
Masaharu Hayashi を著作者とするこの 作品 は クリエイティブ・コモンズの 表示 4.0 国際 ライセンスで提供されています。